library(FFT)
## Loading required package: ggplot2
## Loading required package: Rtsne
## Loading required package: igraph
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
## Loading required package: plotly
##
## Attaching package: 'plotly'
## The following object is masked from 'package:igraph':
##
## groups
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Loading required package: pdfCluster
## pdfCluster 1.0-3
##
## Attaching package: 'pdfCluster'
## The following object is masked from 'package:plotly':
##
## groups
## The following object is masked from 'package:igraph':
##
## groups
## Loading required package: RANN
## Loading required package: flashClust
##
## Attaching package: 'flashClust'
## The following object is masked from 'package:stats':
##
## hclust
## Loading required package: dimRed
## Loading required package: DRR
## Loading required package: kernlab
##
## Attaching package: 'kernlab'
## The following object is masked from 'package:ggplot2':
##
## alpha
## Loading required package: CVST
## Loading required package: Matrix
##
## Attaching package: 'CVST'
## The following object is masked from 'package:igraph':
##
## blocks
##
## Attaching package: 'dimRed'
## The following object is masked from 'package:stats':
##
## embed
## The following object is masked from 'package:base':
##
## as.data.frame
## Loading required package: RColorBrewer
## Loading required package: pheatmap
## Loading required package: colorRamps
## Loading required package: RSpectra
## Loading required package: foreach
## Loading required package: doParallel
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: e1071
## Loading required package: ggraph
## Loading required package: cluster
## Loading required package: umap
## Loading required package: leiden
## conda environment r-reticulate installed
## python modules igraph and leidenalg installed
library(PTA)
## Loading required package: fda
## Loading required package: splines
## Loading required package: fds
## Loading required package: rainbow
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
##
## select
## Loading required package: pcaPP
## Loading required package: RCurl
##
## Attaching package: 'fda'
## The following object is masked from 'package:graphics':
##
## matplot
##
## Attaching package: 'PTA'
## The following object is masked from 'package:FFT':
##
## getMapId
library(leiden)
cluLeiden = NULL
load("../../LISA2new/data/simulation_res.RData")
if(is.null(pcd)){
pcd = pcaL1(data=daC, pvalue=0.05, pcaSel = 1, sel=2, center=TRUE, scale=FALSE, ith=1e-2)
}
if(is.null(Z)){
nc = getPCAN(pcd$pca_out,pvalue=0.05, sel=2, ith=1e-2)
Z = pcd$pca_out$x[,1:nc]
}
dim(Z)
## [1] 1300 10
Compute UMAP using the selected PCs.
library(plotly)
# For plot
set.seed(1)
if(is.null(uY)){
emb2 <- umap::umap(Z)
uY = emb2$layout
}
plotAlls(X=uY,draw=1, t=NULL, cluid=as.character(branch),
type=2,size =6,opacity=1, titlename=paste0("UMAP: PC number = ", ncol(Z)))
ks = seq(4,20, by=2)
if(is.null(cluLeiden)){
tclu3 = proc.time()
cluLeiden = list()
#setLeignPath()
for(i in 1:length(ks)){
cluLeiden[[i]] = cmcluster(Z, ksize = ks[i], cmmethod="leiden")
}
tclu4 = proc.time() - tclu3
}
Plot clustering
lp <- htmltools::tagList()
for (i in 1:length(ks)) {
cluid = cluLeiden[[i]]$cluid
table(cluid)
p1 = plotAlls(X=uY,draw=1, t=NULL, cluid=as.character(cluid), type=2,size =6,opacity=1,
titlename=paste("UMAP,k=", ks[i], sep=""))
lp[[i]] <- as.widget(p1)
}
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
## Warning: 'as.widget' is deprecated.
## Use 'as_widget' instead.
## See help("Deprecated")
lp
k=5
cluid = cluLeiden[[k]]$cluid
table(cluid)
## cluid
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 111 109 104 101 97 97 96 94 92 88 74 65 55 51 36 30
plotAlls(X=uY,draw=1, t=NULL, cluid=as.character(cluid), type=2,size =6,opacity=1,
titlename=paste("UMAP,k=", ks[k], sep=""))
library("netbiov")
plotCommNet(Z, ksize=ks[k], cluid, colname= rownames( brewer.pal.info)[12],
gtype = "Kamada-Kawai", seedN = 1)
centid = getCluCenterID(Z, cluid)
peakid = searchPeaks(gg=cluLeiden[[k]]$knng$g, cluid)
X <- autoLISOMAP(Z, peaks=peakid, ndim = 3, L1=NULL, L2=NULL, N=1)
## 2021-03-16 10:44:21: constructing knn graph
## 2021-03-16 10:44:21: calculating geodesic distances
## 2021-03-16 10:44:21: Classical Scaling
## 2021-03-16 10:44:21: constructing knn graph
## 2021-03-16 10:44:21: calculating geodesic distances
## 2021-03-16 10:44:21: embedding
Visualize 3D L-ISOMAP
PZ = rbind(uY, uY[peakid,])
pc = c(cluid, rep("peaks", length(peakid)))
plotAlls(X=PZ,draw=1, t=NULL, cluid=as.character(pc), type=2,size =6,opacity=1,
titlename=paste("Peaks in UMAP,k=", ks[k], sep=""))
#### 4. Visualization
#plotAlls(X=X,draw=3, t=NULL, cluid=as.character(branch), type=2,size =6,opacity=1, titlename="ISOMAP")
plotAlls(X=X,draw=3, t=NULL, cluid=as.character(cluid), type=2,size =5,opacity=1, titlename="ISOMAP")
##1. Compute neighbor distance of clusters in graph built by PCA
msM = getCluNeighbors(Z, k=ks[k], cluid=cluid)
## 2021-03-16 10:44:22: Get neighbors
## 2021-03-16 10:44:22: End
##2. Compute graph distance matrix
dsM = distGraphClu(Z, ksize=ks[k], cluid=cluid)
## 2021-03-16 10:44:22: Get neighbors
## 2021-03-16 10:44:23: End
rt = 2 # clu
g1 = getMSTLeafRootISO_3(neigDist = msM, cenDist=dsM, cluid=cluid, L2=3, root=rt,
leaves=list(a=c(5,7,"parallel"), b=c(8,9,"parallel"),
c=c(16,12,"linear"), d=10, f=4))
## [1] "L2 is too small to build a connected graph!"
## [1] 4
## [1] 17
clucells = getCluCellTh(cluid=cluid, branch, th=0.2, subN=10)
g2 <- set.vertex.attribute(g1$g5, "name", value=clucells)
ctyM = getCellSM(cluid, branch)
cellN = as.character(names(table(branch)))
g3 <- set.vertex.attribute(g1$g5, "name", value= paste("Y_",1:max(cluid),sep=""))
pieTreeRoot(mst=g3, root=rt, ctyM, cellN, vs=1, ew=1, cex=2, seed=1,colorname="Paired") # Plot the cluster tree
#### Plot 3D trajectory
plot3DTreeCells(X[,1:3], peaks =NULL, cluid= as.character(cluid), celltype = as.character(branch),
colorname="Paired", mst=g1$g5, pointsize = 6,labelsize=15, lineW=10)
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
#plot3DTreeCells(X[,1:3], peaks =NULL, cluid= as.character(cluid),colorname="Paired",
# celltype = as.character(cluid),mst=g1$g5, pointsize = 6,labelsize=15, lineW=10)
#### Plot the tree
Y <- getCenters(X, cluid)
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
cluidS = as.character(cluid)
anno = cbind(cluidS, as.character(branch))
colnames(anno) = c("cluster", "celltype")
anno = as.data.frame(anno)
g3 <- set.vertex.attribute(g1$g5, "name", value= paste("Y_",1:max(cluid),sep=""))
#g3 <- set.vertex.attribute(g1, "name", value= paste("Y_",1:max(cluid),sep=""))
plotCellTree(X=X, uY, cluid=cluid, mst=g3, root=rt, colorby="celltype", anno=anno, height=0.2,
width=0.3, pointSize=2, edgeW=0.5, textSize=6, legendTextS=20, legendColS=10, hv=0.5,
vv=-1.2, colorName = "Paired")
rt = 2# root
g2 <- set.vertex.attribute(g1$g5, "name", value= paste("Y_", 1:max(cluid),sep = ""))
res = simpPTime(Z, cluid=cluid, cenTree=g2, root=rt)
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## [1] "Y_1" "Y_2" "Y_3" "Y_4" "Y_5" "Y_6" "Y_7" "Y_8" "Y_9" "Y_10"
## [11] "Y_11" "Y_12" "Y_13" "Y_14" "Y_15" "Y_16"
# save(a=pcd, b=Z, c=uY, d=cluLeiden, e=X, f=g1, g=branch, h=daC, i=Y, j=res,k=g2, file = "../data/simulation_res.RData")
plotAlls(X=uY,draw=1, t=res$pt, cluid=as.character(branch), type=1,size =6,opacity=1, titlename="UMAP")
#plotAlls(X=X,draw=3, t=res$pt, cluid=as.character(branch), type=1,size =6,opacity=1, titlename="ISOMAP")
User can apply PTA for each lineage 1. User should specify each lineage first as input of PTA 2. To reduce the computation time, user can do average on the cells along the lineage
lpath1 = c(2,15,11,9)
lpath2 = c(2,15,11,8)
lpath3 = c(2,14,13,6,4)
lpath4 = c(2,14,13,6,3,7)
lpath5 = c(2,14,13,6,3,5)
lpath6 = c(2,14,13,1,16,12)
lpath7 = c(2,14,13,1,10)
cluid = cluLeiden[[5]]$cluid
fpath = "../../LISA2new/data/simu_PTA_lineage_1.RData"
if(!file.exists(fpath)){
usePTA(sd=daC, cluid=cluid, path=lpath1, pt = res$pt, fname=fpath,
it = 25, intv = 200, feature.flag = TRUE,
err_flag = 1, scale=FALSE, nCores=1)
}
load(fpath)
library(reshape2)
library(colorRamps)
library(grid)
library(gridExtra)
library(pheatmap)
singlePath = c(2,15,11,9)
cid = singlePath
da = NULL
for(i in 1:length(cid)){
k= which(cluid==cid[i])
pt = res$pt[k]
od = order(pt)
k = k[od]
da = cbind(da, daC[,k])
}
plotHeatMap <- function(da1,tv, fz=20, rowsize = 6, colsize=4,barsize=10, showRname=F, showCname=F,
colname=NULL, title=""){#
#cs = colorRampPalette(rev(c("#D73027", "#FC8D59", "#FEE090",
# "#FFFFBF", "#E0F3F8", "#91BFDB", "#4575B4")))(100);
if(is.null(colname)){
#cs = colorRampPalette(rev(c("#D73027", "#FC8D59", "#FEE090",
# "#FFFFBF", "#E0F3F8", "#91BFDB", "#4575B4")))(100);
cs = colorRampPalette(c("#0392cf","#7bc043", "#fdf498","#ffff66","#ff3300"))(50)
}else{#colname="Spectral"
cs = colorRampPalette(brewer.pal(brewer.pal.info[colname,1],colname))(100)
}
a = pheatmap(da1, color = cs, breaks = NA, scale = "none", cluster_rows = F,
cluster_cols = F, legend = TRUE, annotation = NA,
annotation_colors = NA, annotation_legend = TRUE, filename = NA, width = NA,
height = NA, main = title,
show_rownames=showRname,show_colnames=showCname,
silent=TRUE, fontsize = fz,fontsize_row = rowsize, fontsize_col = colsize)
if(tv[1] < 0){
od = order(tv, decreasing = TRUE)
tv = tv[od]
}else{
od = order(tv)
tv = tv[od]
}
df <- data.frame(top_v=1:length(tv), loadings=tv)
p<- ggplot(data=df, aes(x=top_v, y=loadings)) +
geom_bar(stat="identity") + coord_flip() +
theme(panel.background = element_blank(), axis.text=element_text(size=barsize),
axis.title=element_text(size=barsize,face="bold"),
axis.ticks.x=element_blank())+
scale_x_continuous(position = "top") +
scale_y_continuous(position = "right")
#+ theme_classic()
#+ theme_minimal()
lm=rbind(c(1,1,1,1,1,2))
grid.arrange(grobs = list(a[[4]],p), layout_matrix = lm)
}
str="Main trends in simulation lineage 2->9"
feature.flag=TRUE
titlefontSize = 2
# out$uName
##################all trends
mout = out
timevec = 1:nrow(mout$B)
S1 = mout$B %*% mout$beta
f1 <- splinefun(timevec, S1)
T = length(timevec)
max1 = max(mout$prd)
min1 = min(mout$prd)
mout = out$rank2
S2 = mout$B %*% mout$beta
f2 <- splinefun(timevec, S2)
max2 = max(mout$prd)
min2 = min(mout$prd)
mout = out$rank3
S3 = mout$B %*% mout$beta
f3 <- splinefun(timevec, S3)
max3 = max(mout$prd)
min3 = min(mout$prd)
ymin = min(c(S1,S2,S3))
ymax = max(c(S1,S2,S3))
vmax = max(max1, max2, max3)
vmin = min(min1, min2, min3)
plot(x= timevec, y= S1, ylim=c(ymin, ymax), type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black", main= str, cex.main= titlefontSize)
curve(f1(x), from=timevec[1], to=timevec[T], col="red", lty=1,lwd=2, add=TRUE)
points(x= timevec, y= S2, type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black")
curve(f2(x), from=timevec[1], to=timevec[T], col="blue", lty=1,lwd=2, add=TRUE)
points(x= timevec, y= S3, type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black")
curve(f3(x), from=timevec[1], to=timevec[T], col="green", lty=1,lwd=2, add=TRUE)
par(xpd = T)
legend(1, ymax+1, title = " ", legend=c('Rank1','Rank2','Rank3'),
col = c('red','blue','green'), horiz = F,bty='n',
pch = 14,pt.cex = 2, cex=1.5,ncol= 10 )
plot(x= timevec, y= S1, ylim=c(ymin, ymax), type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black", main= str, cex.main= titlefontSize)
curve(f1(x), from=timevec[1], to=timevec[T], col="red", lty=1,lwd=2, add=TRUE)
points(x= timevec, y= S2, type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black")
curve(f2(x), from=timevec[1], to=timevec[T], col="blue", lty=1,lwd=2, add=TRUE)
par(xpd = T)
legend(1, ymax+4, title = " ", legend=c('Rank1','Rank2' ),
col = c('red','blue'), horiz = F,bty='n',
pch = 16,pt.cex = 2, cex=1.5,ncol= 10 )
PTA.plotCV(err,feature.flag = feature.flag)
plot(mout$a[,1],type="p", xlab="genes", ylab="a", lwd=2, font=2)
TH1 = 0.03
corTH = 0
mout = out
na <- c(length(mout$a[mout$a>TH1]), length(mout$a[mout$a<-TH1]),length(mout$a[abs(mout$a)<TH1]) )
names(na) = c(paste('a>',TH1),
paste('a<',-TH1),
paste('|a|<',TH1) )
na
## a> 0.03 a< -0.03 |a|< 0.03
## 39 0 0
gid = 1:length(mout$a)
uid = gid[mout$a>TH1]
did = gid[mout$a<-TH1]
# up and down expressed marker genes
upgene1 = as.character(out$uName[uid])
downgene1 = as.character(out$uName[did])
############################ Plot Trend
timevec = 1:nrow(mout$B)
S = mout$B %*% mout$beta
f1 <- splinefun(timevec, S)
T = length(timevec)
plot(x= timevec, y= S, type="p", xlab="Time", ylab="S", lwd=2, font=2, col="black",ylim=range(S), main="Trend")
curve(f1(x), from=timevec[1], to=timevec[T], col="red", lty=1,lwd=2, add=TRUE)
################################ predicted marker
mout = out
rownames(mout$prd) = mout$uName
nid = which( abs(mout$a[,1]) > TH1)
uWt = mout$a[nid,1]
od = order(uWt,decreasing = TRUE)
uid = nid[od]
da1 = mout$prd[uid,]
plotHeatMap(da1, uWt[od], fz=20, title = "Predicted data gene rank 1")
x = mout$x[,1,]
rownames(x) = mout$uName
nid = which( abs(mout$a[,1]) > TH1)
uWt = mout$a[nid,1]
od = order(uWt,decreasing = TRUE)
uid = nid[od]
da1 = x[uid,]
cuWt = uWt[od]
corS= apply(da1, 1, function(x) {cor(x, S1) })
id = which(abs(corS)>=corTH)
plotHeatMap(da1[id,], cuWt[id], fz=20, title = "Raw data gene rank 1")
kid = getMapId(rownames(da), rownames(da1[id,]))
plotHeatMap(da[kid,], cuWt[id], fz=20, showRname=T, title = "Raw data gene rank 1")
TH1 = 0.03
mout = out$rank2
############################ plot loadings
PTA.plotCV(out$err2,feature.flag = feature.flag)
plot(mout$a[,1],type="p", xlab="genes", ylab="a", lwd=2, font=titlefontSize)
na <- c(length(mout$a[mout$a>TH1]), length(mout$a[mout$a<-TH1]),length(mout$a[abs(mout$a)<TH1]) )
names(na) = c(paste('a>',TH1),
paste('a<',-TH1),
paste('|a|<',TH1) )
na
## a> 0.03 a< -0.03 |a|< 0.03
## 23 0 0
gid = 1:length(mout$a)
uid = gid[mout$a>TH1]
did = gid[mout$a<-TH1]
# up and down expressed marker genes
upgene2 = as.character(out$uName[uid])
downgene2 = as.character(out$uName[did])
############################ Plot Trend
timevec = 1:nrow(mout$B)
S = mout$B %*% mout$beta
f1 <- splinefun(timevec, S)
T = length(timevec)
plot(x= timevec, y= S, type="p", xlab="Time", ylab="S", lwd=2, font=titlefontSize, col="black",ylim=range(S), main="Trend")
curve(f1(x), from=timevec[1], to=timevec[T], col="red", lty=1,lwd=2, add=TRUE)
Plot
################################ predicted marker
mout = out$rank2
rownames(mout$prd) = out$uName
nid = which( abs(mout$a[,1]) >TH1)
uWt = mout$a[nid,1]
od = order(uWt,decreasing = TRUE)
uid = nid[od]
da1 = mout$prd[uid,]
plotHeatMap(da1, uWt[od], fz=20, title = "Predicted data gene rank 2")
x = out$x[,1,]
rownames(x) = out$uName
nid = which( abs(mout$a[,1]) >TH1)
uWt = mout$a[nid,1]
od = order(uWt,decreasing = TRUE)
uid = nid[od]
da1 = x[uid,]
corS= apply(da1, 1, function(x) {cor(x, S2) })
id = which(abs(corS)>=corTH)
cuWt = uWt[od]
plotHeatMap(da1[id,], cuWt[id], fz=20, title = "Raw data gene rank 2")
kid = getMapId(rownames(da), rownames(da1[id,]))
plotHeatMap(da[kid,], cuWt[id], fz=20, showRname=T, title = "Raw data gene rank 2")